function sac3()

QueMessage('SAC3 Starting', 1);

sfm = getmainselection;

if(sfm > 0)
    pflag = getplotflag;
    QueMessage('SAC3 analysis, Sliding windows', 1); % clear the que
    for i = 1:length(sfm)
        sac3_go(sfm(i), pflag);
    end;
end;

function sac3_go(sf, pflag)

global DFILE ALLCH VOLTAGE CONTROL



if(pflag && ~isempty(CONTROL(sf).SAC3))
    print_sac3(sf, pflag);
    return;
end;

isrc = 1;
t0 = 0;
tmax = 1000;
spike_thresh = 0;

X=[];
if(~isempty(CONTROL(sf).spike))
    X=getstructs(sf, 'spike.latency');
    X=X{:};
    ns = 0;
    for i = 1:length(X)
        ns = ns + length(X{i});
    end;
    fprintf(1, 'spike data for %s contains: %d trials with %d spikes\n', ...
        CONTROL(sf).filename, length(X), ns);
end;
if(isempty(X) | isempty(CONTROL(sf).spike)) % only read spikes if we really need to - otherwise the data is in the database already
    fprintf(1, 'Spike data not available: Reading %s\n', CONTROL(sf).filename);
    [DFILE, DPAR, err] = ana_setup(DFILE, sf);
    if(err ~= 0)
        return;
    end;

    if(~isempty(ALLCH))
        VOLTAGE = ALLCH{isrc};
    end;
    %[protocol, rate, records, pts, frec, lrec, time, TM, ZT, TL, VL] = analysis_setup2(DFILE, sf);
    [first_spike, first_isi, nr_spikes, splat]=find_spikes2(DFILE, VOLTAGE, t0, tmax, spike_thresh);
    X = {splat.latency};
    spike.latency = {splat.latency};
    spike.source = {splat.source};
    spike.fsl = first_spike;
    CONTROL(sf).spike=spike;
end;


binw = 0.025;
nwin = 10;
winlen = 100; % msec
winstep = 100; % step (if less than winlen, we are using a "sliding" window)
twin = winlen;
for n = 1:nwin
    start = (n - 1) * winstep;
    dur = winlen;
    tmid(n) = start + dur/2;
    [y, yh1, hx1, mr1] = sac(X, twin, binw, start, dur);
    % perform gaussian fits on the histograms to identify the peaks and measure
    % them. Using Molitor's routines (mrqfit).
    % [...] = MRQFIT(F, P, X, Y, SIG, VP, LB, UB, IMAX, TOL)
    % F = 'gaussian', P = [A0 A1 M1 S1 A2 M2 S2 ... AN MN SN]
    % SIG is sigma Y (def = 1); VP is vary array per parameter [0 or 1];
    % LB, UB are upper and lower bounds.
    %
    if(isempty(y))
        hx1c{n} = [];
        fp11 = NaN*ones(4, 1);
        chisq11 = NaN;
        fitc11 = fp11;
        err11 = fp11;
        dep11 = fp11;
        niter11 = 0;
        fwhm = NaN;
        ci = NaN;
        hx1 = [];
        yg1 = [];
        yh1 = [];
        mr1 = [];
        fprintf(1, 'Sac3 (%d) y is empty\n', n);
    else

        hx1c{n} = hx1+0.5*(hx1(2)-hx1(1));
        gpar1 = [1, 5, 0, 1]; % single gaussian centered on 0
        scf = 1;
        gulim1 = [0, 100, 10*scf, 20*scf];
        gllim1 = [0, 0, 0.0, 0.01];
        gvar1 =  [0 1 0 1];
        nitermax = 100;

        [fp11, chisq11, niter11, fitc11, err11, dep11] = mrqfit('gaussian', gpar1, hx1c{n}*scf, yh1, [], gvar1, gllim1, gulim1, nitermax, []);



        yg1 = gaussfunc(hx1c{n}, fp11);
        fwhmfac = 2*sqrt(2*log(2)); % note - log is ln.
        fwhm = fwhmfac*fp11(4); % full width at half maximal height.
        ci = max(yg1);
        fprintf(1, 'Sac3 (%d) t = %.1f  ci = %.2f\n', n, tmid(n), ci);
    end;

    % generate the results structure: SAC
    SAC3.hx1{n} = hx1; % save the histograms.
    SAC3.hy1{n} = yh1;
    SAC3.yg1{n} = yg1;
    SAC3.mr1{n} = mr1;

    % gaussian fit results:
    SAC3.Gfit1{n} = fp11; % parameters
    SAC3.niter1(n) = niter11; % iterations
    SAC3.chisq1(n) = chisq11; % fit error
    SAC3.err1{n} = err11; % parameter estimate error
    SAC3.dep1{n} = dep11; % dependency between parameters
    SAC3.fwhm1(n) = fwhm;
    SAC3.NPH1(n) = ci;
end;
SAC3.start = start;
SAC3.dur = dur;
SAC3.nwin = nwin;
SAC3.tmid = tmid;
SAC3.mr1 = mr1;
SAC3.ns = ns;
SAC3.twin = twin;

%

CONTROL(sf).SAC3 = SAC3; % save in the database.

print_sac3(sf, pflag);


function print_sac3(sf, pflag)
%
%
% function just to print the current SAC -
%

global CONTROL;
SAC3 = CONTROL(sf).SAC3;

hf = findobj('tag', 'SAC3');
if(isempty(hf))
    hf = figure;
    set(hf, 'Tag', 'SAC3');
end;



figure(hf);
clf;
set(hf, 'color', [1 1 1]); % white background...

% text area
subplot('position', [0.1, 0.90, 0.8, 0.095]);
axis([0,1,0,1])
axis('off')
ht(1)=text(0,0.80,sprintf('%-12s R[%d:%d]     %-8s  [%s]', ...
    CONTROL(sf).filename, CONTROL(sf).recbeg, CONTROL(sf).recend, ...
    CONTROL(sf).protocol, date), 'Fontsize', 10);
set(ht(1), 'interpreter', 'none'); % un-TeX the line - this is a filename and won't have tex chars, but might have an underscore.
text(0,0.6,sprintf('Solution:%-12s  gain:%4.1f ', CONTROL(sf).solution, CONTROL(sf).igain), 'FontSize', 8);
text(0,0.4,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
    CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 8);
text(0,0.200,sprintf('Windows: %.1f-%.1f (nwin=%d)  mean rate: %.2f s/s tot spikes: %d', ...
    SAC3.start, SAC3.dur, SAC3.nwin, SAC3.mr1, SAC3.ns), 'FontSize', 8);
%text(0, 0.000, sprintf('G1: A=%.2f (%.2f) S=%.3f (%.3f)  NPH = %.2f  FWHM1 = %.3f', ...
%    fp11(2), err11(2), fp11(4), err11(4), SAC.NPH1, SAC.fwhm1), 'Fontsize', 8);

subplot('position', [0.1, 0.075, 0.35, 0.320]);
%bar(sqrt(hx1), yh1, 'histc');
for n = 1:SAC3.nwin
    hx1 = SAC3.hx1{n};
    if(~isempty(hx1))
        hx1c = hx1+0.5*(hx1(2)-hx1(1));
        hy1=SAC3.hy1{n};
        semilogx(hx1c, hy1, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
        hold on
        semilogx(hx1c, SAC3.yg1{n}, 'r');
    end;
    %semilogx(hx1c, fitc11, 'g');
end;
u1 = get(gca, 'Ylim');
h1 = gca;
set(gca, 'Xlim', [0 SAC3.twin]);
xlabel('Delay (ms)');
ylabel('Correlation Index');

if(~isempty(CONTROL(sf).spike))
    X=getstructs(sf, 'spike.latency');
    X=X{:};
    subplot('position', [0.55, 0.075, 0.35, 0.320]); % suplot for RASTER...

for i = 1:length(X)
    plot(X{i}, i*ones(length(X{i}), 1), 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
    if(i > 1) hold on;
    end;
end;
xlabel('Time(ms)');
ylabel('Trial #');
end;

subplot('position', [0.55, 0.480, 0.35, 0.320]);
plot(SAC3.tmid, SAC3.NPH1, 'ks-', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
xlabel('Delay (ms)');
ylabel('Correlation Index');

subplot('position', [0.10, 0.480, 0.35, 0.320]);
plot(SAC3.tmid, SAC3.fwhm1, 'bo-', 'MarkerFaceColor', 'b', 'MarkerSize', 3.5);
xlabel('Delay (ms)');
ylabel('FWHM (ms)');

figure(hf); % bring to the front

if(pflag)
    orient landscape
    [p f] = fileparts(CONTROL(sf).filename);
    pfile = slash4OS(sprintf('%s/Notebook/Sac3/%s-%s-%d', pwd, f, CONTROL(sf).E_C, min(seqparse(CONTROL(sf).reclist))));
    print(hf, '-dpng', pfile);
end;


function [y] = gaussfunc(x, fp)
%
% calculate a gaussian based on x and FP
%
y = fp(1) + (fp(2)/(fp(4)*sqrt(2*pi)))*exp(-((x-fp(3)).^2)/(2*fp(4)^2));

